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Abstract 

We give algorithms for the computation of the d - th digit of certain transcendental 
numbers in various bases. These algorithms can be easily implemented (multiple 
precision arithmetic is not needed), require virtually no memory, and feature run 
times that scale nearly linearly with the order of the digit desired. They make it fea- 
sible to compute, for example, the billionth binary digit of log(2) or % on a modest 
workstation in a few hours run time. 

We demonstrate this technique by computing the ten billionth hexadecimal digit of 
7t, the billionth hexadecimal digits of 7t 2 , log(2) and log 2 (2), and the ten billionth 
decimal digit of log(9/10). 

These calculations rest on the observation that very special types of identities exist 
for certain numbers like 7t, 7t 2 , log(2) and log 2 (2). These are essentially polyloga- 
rithmic ladders in an integer base. A number of these identities that we derive in this 
work appear to be new, for example the critical identity for 7t: 

Y i r 4 2 i i_J\ 

n “ ^I6'l8i+l 8 j + 4 8i + 5 81 + 6 J. 

r = 0 
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1. Introduction. 


It is widely believed that computing just the d-th digit of a number like 7 r is really 
no easier than computing all of the first d digits. From a bit complexity point of 
view this may well be true, although it is probably very hard to prove. What we will 
show is that it is possible to compute just the d-th digit of many transcendentals in 
(essentially) linear time and logarithmic space. So while this is not of fundamentally 
lower complexity than the best known algorithms (for say n or log(2)), this makes 
such calculations feasible on modest workstations without needing to implement 
arbitrary precision arithmetic. 

We illustrate this by computing the ten billionth hexadecimal digit of 7r, the billionth 
hexadecimal digits of tt 2 , log(2) and log 2 (2) , and the ten billionth decimal digit of 
log(9/10). Details are given in Section 4. A previous result in this same spirit 
is the Rabinowitz-Wagon “spigot” algorithm for ir. In that scheme, however, the 
computation of the digit at position n depends on all digits preceding position n. 

We are interested in computing in polynomially logarithmic space and polynomial 
time. This class is usually denoted SC (space = log°^(d) and time = d 0 ^ 1 ) where 
d is the place of the “digit” to be computed). Actually we are most interested in 
the space we will denote by SC* of polynomially logarithmic space and (almost) 
linear time (here we want the time = 0(d log°^(d))). There is always a possible 
ambiguity when computing a digit string base b in distinguishing a sequence of 
digits a(b — 1)(6 — 1)(6 — 1) from (a 4- 1)000. In this particular case we consider 
either representation as an acceptable computation. In practice this problem does 
not arise. 

It is not known whether division is possible in SC, similarly it is not known whether 
base change is possible in SC. The situation is even worse in SC*, where it is not 
even known whether multiplication is possible. If two numbers are in SC* (in the 
same base) then their product computes in time = 0(d 2 log°^ 1 \d)) and is in SC 
but not obviously in SC*. The d 2 factor here is present because the logarithmic 
space requirement precludes the usage of advanced multiplication techniques, such 
as those based on FFTs. 

We will not dwell on complexity issues except to point out that different algorithms 
are needed for different bases (at least given our current ignorance about base 
change) and very little closure exists on the class of numbers with d-th digit com- 
putable in SC. Various of the complexity related issues are discussed in [6,8,9,11,14]. 

As we will show in Section 3, the class of numbers we can compute in SC* in base 
b includes all numbers of the form 


(i.i) 


4 pW 

w 


where p and q are polynomials w r ith integer coefficients and c is a positive integer. 
Since addition is possible in SC*, integer linear combinations of such numbers are 
also feasible (provided the base is fixed). 
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The algorithm for the binary digits of 7T, which also shows that 7T is in SC* in base 
2, rests on the following remarkable identity: 

Theorem 1. The following identity holds: 


00 1 4 

(L2 > ' = El^(87TT-8i 

1=0 


1 


1 


)• 


: + 4 8i + 5 81 + 6 

This can also be written as: 

00 

(1.3) * = EdriF< [p>] = [4, 0, 0, —2, —i, —i,o, oj 


1 = 1 


where the overbar notation indicates that the sequence is periodic . 


Proof This identity is equivalent to: 

/ 1/ ' /5 4 - 8x 3 - 4\/2x 4 - 8a: 5 , 

( 1 4 ) ^ 


which on substituting y := y/2x becomes 




16 y — 16 


2 y 3 + 4 y - 4 


dy. 


The equivalence of (1.2) and (1.4) is straightforward. It follows from the identity 


r i/%/2 r k-\ fi/v* ~ 

/ - -dx= I Yx k ~ 1+ *'dx 

Jo i -* 8 Jo frt 

l 00 l 

= ^S i6 ‘( 8i+fc ) 


r 1 /n/2 


That the integral (1.4) evaluates to n is an exercise in partial fractions most easily 
done in Maple or Mathematica. □ 


This proof entirely conceals the route to discovery. We found the identity (1.2) by 
a combination of inspired guessing and extensive searching using the PSLQ integer 
relation algorithm [3,12]. 

Shortly after the authors originally announced the result (1.2), several colleagues, 
including Helaman Ferguson, Tom Hales, Victor Adamchik, Stan Wagon, Donald 
Knuth and Robert Harley, pointed out to us other formulas for 7r of this type. One 
intriguing example is 

^ 1 , 2 2 1 1/2 1/2 l/4_ \ 

Y I6« ^8i + 1 + 4z + 2 + 4i + 3 4i + 5 4i + 6 4i + 7'’ 
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which can be written more compactly as 




A 1 \ A ■ ‘ - ' 


1 


t=0 


4i -f 1 4i -f 2 4i + 3 


In [2], this and some related identities are derived using Mathematica. 

As it turns out, these other formulas for i r can all be written as formula (1.2) plus 
a rational multiple of the identity 


o = T—( ~ 8 


+ 


+ 


+ 


8 


2 


1 


i = 0 


16* 4* 1 8i + 2 8i + 3 8i + 4 8z‘ + 5 8i + 6 8* + 


7>' 


The proof of this identity is similar to that of Theorem 1. 

The identities of the next section and Section 5 show that, in base 2, 7 r 2 , log 2 (2) 
and various other constants, including {log(2), log(3), . . . ,Iog(22)} are in SC*. (We 
don’t know however if log(23) is even in SC.) 

We will describe the algorithm in the Section 3. Complexity issues are discussed 
in [3,5,6,7,8,9,14,19,21] and algorithmic issues in [5,6,7,8,14]. The requisite special 
function theory may be found in [1,5,15,16,17,20]. 


2. Identities. 

As usual, we define the m-th poly logarithm L m by 

OO 

( 2 -i) ^=y=Y.h' m< *• 

i— 1 


The most basic identity is 

(2.2) — log(l — 2 -n ) = Li(l/2") 

which shows that log(l — 2” n ) is in SC* base 2 for integer n. (See also section 5.) 
Much less obvious are the identities 

(2.3) 7T 2 = 36L 2 (l/2) - 36L 2 (l/4) - 12L 2 (l/8) + 6L 2 (l/64) 

and 

(2.4) l°g 2 (2) = 4L 2 (l/2) - 6L 2 (l/4) - 2L 2 (l/8) + L 2 (l/64) . 

These can be written as 

OO 

(2.5) *2 = 36^-^, [a,] = [1,-3, -2, -3, 1,0] 

1 = 1 
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(2.6). 


[bi] = [2, -10, -7, -10, 2, -lj . 


°o , 

I !o S ! (2) = 2£^, 

1=1 

Here the overline notation indicates that the sequences repeat. Thus we see that 
7 r 2 and log 2 (2) are in SC* in base 2. These two formulas can alternately be written 

.16 24 8 6 1 , 

(6i + l) 2 ~ (6i + 2) 2 - (6i + 3) 2 (6f + 4) 2 + (6i + 5) 2 ' 

2/ox IV- x / -16 16 40 14 10 1 . 

° S 2 ~ 8 h> 64< ( 6 0 2 + (6*+ l) 2 (6i + 2) 2 (6i + 3) 2 (6i + 4) 2 + (6f + 5) 2 >' 

Identities (2.3)-(2-6) are examples of polylogarithmic ladders in the base 1/2 in 
the sense of [16]. As with (1.2) we found them by searching for identities of this 
type using an integer relation algorithm. We have not found them directly in print. 
However (2.5) follows from equation (4.70) of [15] with a = tt/3,/? = 7t/ 2 and 7 = 
7t/ 3. Identity (2.6) now follows from the well known identity 

(2.7) 12L 2 (l/2) = 7r 2 -61og 2 (2). 

A distinct but similar formula that we have found for 7r 2 is 


which can be derived from the methods of section 1. 

There are several ladder identities involving L3: 

(2.8) 35/2C(3) - 7T 2 log(2) = 36L 3 (l/2) - 18Z, 3 (l/4) - 4L 3 (l/8) + A 3 (l/64) , 

(2.9) 21og 3 (2) - 7C(3) = — 24L 3 (l/2) + 18L 3 (l/4) + 4L 3 (l/8) - I 3 (l/64) , 

(2.10) 101og 3 (2)-27r 2 log(2) = -48L 3 (l/2) + 54L 3 (l/4) + 12L 3 (l/8)-3L 3 (l/64). 


( 


16 


16 


8 


16 


(8i + l) 2 (Si + 2) 2 (8i + 3) 2 (8i + 4) 2 (8i + 5) 2 (Si + 6) 2 (8i + 7) : 


■)- 


00 1 


TT 2 - 9 V 1 

8 2^ 64' 

j = 0 


The favored algorithms for 7r of the last centuries involved some variant of Machines 
1706 formula: 

. 7T 1 1 

(2.11) — = 4 arctan - — arctan . 

4 0 Zoa 


There are many related formula [15,16,17,20] but to be useful to us all the arguments 
of the arctans have to be a power of a common base, and we have not discovered 
any such formula for 7 r . One can however write 


( 2 . 12 ) 


— = 2 arctan — = -f arctan — = 

2 y/2 Vs 
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This can be written as 


(2.13) V2tt = 4/(1/2) + /(1/8) where /(*) := £ 

2=1 


and allows for the calculation of \pi it in SC*. 

Another two identities involving Catalan’s constant G, it and log(2) are: 


(2.14) 

and 

(2.15) 


G - 


w l°g(2) 


= E 


i~ 1 




5 2 log (2) 

96 8 


t=l 


[Ci] = [1, 1, 1,0, -1. -1,-1, OJ 


°° d 

= E^nr7^' [ 4 ] = [i,o, -i,-i, -i,o, i,i] 

/l 2 -U- 


These may be found in [17 p. 105, p. 151]. Thus SG — 7rlog(2) is also in SC* in 
base 2, but it is open and interesting as to whether G is itself in SC* in base 2. 

A family of base 2 ladder identities exist: 


Ml/64) M 1/8) 2 L m (l/4) 4 L m (l/2) 5(-log(2)) m 

6 m_l 3 m ^ 1 2 m_1 9 9m! 


7r 2 (-log(2)) m - 2 ^(-log(2)r- 4 403C(5)(— log(2)) m - 5 

54 (m - 2)! 486 (m - 4)! 1296 (m - 5)! 

The above identity holds for 1 < rn < 5; when the arguments to factorials are 
negative they are taken to be infinite so the corresponding terms disappear. See 
[16, p. 45], 


As in the case of formula (L2) for 7 r, colleagues of the authors have subsequently 
pointed out several other formulas of this type for various constants. Three exam- 
ples reported by Knuth, which are based on formulas in [13, p. 17, 18, 22, 47. 139], 
are 


^'"<l + ^) = Ei(i7TT + 


1 


1 = 0 


16* v 8z + 1 


V2 1/4 1/8 

8* 4“ 3 Si + 5 82 — 7 


00 1 1 

V / 2arctan(l/v / 2) = ^ 


2=0 


16* K Si + 1 


1/2 1/4 1/8 . 

82 -f- 3 82 4* 5 82 ~b 7 


1 

arctan(l/3) = ( 

i=o 1 


1 

8 ?' + 1 


1 1/2 1/4 . 

8i + 2 8?' + 4 82 -f 5 


Thus these constants are also in class SC*. Some other examples can be found in 

[ 18 ]. 
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3. The Algorithm. 

Our algorithm to compute individual base-6 digits of certain constants is based on 
the binary scheme for exponentiation, wherein one evaluates x n rapidly by succes- 
sive squaring and multiplication. This reduces the number of multiplications to 
less than 21og 2 (n). According to Knuth [14], where details are given, this trick 
goes back at least to 200 B.C. In our application, we need to perform exponentia- 
tion modulo a positive integer c, but the overall scheme is the same — one merely 
performs all operations modulo c. An efficient formulation of this algorithm is as 
follows. 

To compute r = b n mod c, first set t to be the largest power of two < n, and set 
r = 1. Then 


A: if n > t then r br mod c; n «— n — t\ endif 
t <- i/2 

if t > 1 then r <— r 2 mod c; go to A; endif 

Here and in what follows, “mod” is used in the binary operator sense, namely as the 
binary function defined by x mod y := x — [ x/y)y . Note that the above algorithm 
is entirely performed with positive integers that do not exceed c 2 in size. Thus it 
can be correctly performed, without round-off error, provided a numeric precision 
of at least 1 + 2 log 2 c bits is used. 

Consider now a constant defined by a series of the form 



1 

b ck p(k) ’ 


where 6 and c are positive integers and p(k) is a polynomial with integer coefficients. 
First observe that the digits in the base 6 expansion of ^beginning at position n + 1 
can be obtained from the fractional part of b n S-. Thus we can write 

(3. 


4) 


°° Ln — ck 

b n S~ mod 1 = V — — — mod 1 
to ><*> 


= E 


k=0 


b n ck mod p(fc) 
p(k) 


mod 1 + 


E 


fc = Ln/cJ + l 


b n-ck 

p(k) 


mod 1 


For each term of the first summation, the binary exponentiation scheme is used 
to evaluate the numerator. Then floating-point arithmetic is used to perform the 
division and add the result to the sum mod 1. The second summation, where 
the exponent of 6 is negative, may be evaluated as written using floating-point 
arithmetic. It is only necessary to compute a few terms of this second summation, 
just enough to insure that the remaining terms sum to less than the “epsilon” of 
the floating-point arithmetic being used. The final result, a fraction between 0 and 
1, is then converted to the desired base 6. 
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Since floating-point arithmetic is used here in divisions and in addition modulo 1, 
the result is of course subject to round-off error. If the floating-point arithmetic 
system being used has the property that the result of each individual floating-point 
operation is in error by at most one bit (as in systems implementing the IEEE 
arithmetic standard), then no more than log 2 (2n) bits of the final result will be 
corrupted. This is actually a generous estimate, since it does not assume any 
cancelation of errors, which would yield a lower estimate. In any event, it is clear 
that ordinary IEEE 64-bit arithmetic is sufficient to obtain a numerically significant 
result for even a large computation, and “quad precision” (i.e. 128-bit) arithmetic, 
if available, can insure that the final result is accurate to several digits beyond 
the one desired. One can check the significance of a computed result beginning 
at position n by also performing a computation at position n + 1 or n — 1 and 
comparing the trailing digits produced. 

The most basic interesting constant whose digits can be computed using this scheme 
is 

00 1 

b s< 2 > = E w 

k — 1 

in base 2. Using this scheme to compute hexademical digits of 7 r from identity (1.2) 
is only marginally more complicated, since one can rewrite formula (1.2) using four 
sums of the required form. Details are given in the next section. In both cases, in 
order to compute the n-th binary digit (or a fixed number of binary digits at the 
n-th place) we must sum 0(n ) terms of the series. Each term requires 0(log(n)) 
arithmetic operations and the required precision is O(log(n)) digits. This gives 
a total bit complexity of 0(n log(n)A/(log(n))) where M(j) is the complexity of 
multiplying j bit integers. So even with ordinary multiplication the bit complexity 
is 0(n log 3 (n)). 

This algorithm is, by a factor of log(log(log(n))), asymptotically slower than the 
fastest known algorithms for generating the n-th digit by generating all of the first 
n digits of log(2) or 7r [7]. The asymptotically fastest algorithms for all the first 
n digits known requires a Strassen-Schonhage multiplication [19]; the algorithms 
actually employed use an FFT based multiplication and are marginally slower than 
our algorithm, from a complexity point of view, for computing just the n-th digit. Of 
course this complexity analysis is totally misleading: the strength of our algorithm 
rests mostly on its easy implementation in standard precision without requiring 
FFT methods to accelerate the computation. 

It is clear that the above methods can easily be extended to evaluate digits of 
contstants defined by a formula of the form 



¥ k q(k)' 


where p and q are polynomials with integer coefficients and c is a positive integer. 
Similarly if p and q are slowly growing analytic functions of various types the 
method extends. 
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4. Computations. 

We report here computations of w, log(2), log 2 (2), tt 2 and log(9/10), based on 
the formulas (1.1), (2.2), (2.5), (2.6) and the identity log(9/10) = -Li(l/10), 
respectively. 

Each of our computations employed quad precision floating-point arithmetic for 
division and sum mod 1 operations. Quad precision is supported from Fortran on 
the IBM RS6000/590 and the SGI Power Challenge (R8000), which were employed 
by the authors in these computations. We were able to avoid the usage of explicit 
quad precision in the exponentiation scheme by exploiting a hardware feature com- 
mon to these two systems, namely the 106-bit internal registers in the multiply-add 
operation. This saved considerable time, because quad precision operations are 
significantly more expensive than 64-bit operations. 

Computation of 7 r 2 and log 2 (2) presented a special challenge, because one must 
perform the exponentiation algorithm modulo k 2 instead of k . When n is larger than 
only 2 13 , some terms of the series (2.5) and (2.6) must be computed with a modulus 
k 2 that is greater than 2 26 . Squares that appear in the exponentiation algorithm will 
then exceed 2 52 , which is the nearly the maximum precision of IEEE 64-bit floating- 
point numbers. When n is larger than 2 2 ®, then squares in the exponentiation 
algorithm will exceed 2 104 , which is nearly the limit of quad precision. 

This difficulty can be remedied using a method which has been employed for ex- 
ample in searches for Wieferich primes [10]. Represent the running value r in the 
exponentiation algorithm by the ordered pair (r*i,r 2 ), where r = r x + kr 2 , and 
where r x and r 2 are positive integers less than k. Then one can write 

r 2 = (r x + fcr 2 ) 2 = r 2 + 2r x r 2 fc + r\k 2 

When this is reduced mod k 2 , the last term disappears. The remaining expression is 
of the required ordered pair form, provided that r\ is first reduced mod k, the carry 
from this reduction is added to 2r x r 2 , and this sum is also reduced mod k. Note 
that this scheme can be implemented with integers of size not exceeding 2k 2 . Since 
the computation of r 2 mod k 2 is the key operation of the binary exponentiation 
algorithm, this means that ordinary IEEE 64-bit floating-point arithmetic can be 
used to compute the n-th hexadecimal digit of 7r 2 or log 2 (2) for n up to about 
2 24 . For larger n, we still used this basic scheme, but we employed the multiply- 
add “trick” mentioned above to avoid the need for explicit quad precision in this 
section of code. 

Our results are given below. The first entry, for example, gives the 10 6 -th through 
10 6 -(- 13-th hexadecimal digits of 7r after the “decimal” point. In all cases we did the 
calculations twice — the second calculation was similar to the first, except shifted 
back one position. Since this changes all the arithmetic performed, it is a highly 
rigorous validity check. Thus we believe that all the digits shown below are correct. 
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Constant: 

Base: 

Position: 

Digits from Position: 

7 r 

16 

10 6 

26C65E52CB4593 



10 7 

17AF5863EFED8D 



10 8 

ECB840E21926EC 



10 9 

85895585 A0428B 



10 10 

921C73C6838FB2 

I°g(2) 

16 

10 6 

418489A9406EC9 



10 7 

815F479E2B9102 



10 s 

E648F40940E13E 



10 9 

B1EEF1252297EC 

7T 2 

16 

10 6 

685554E1 228505 



10 7 

9862837AD8AABF 



10 8 

4861AAF8F861BE 



10 9 

437A2BA4A13591 

l°g 2 (2) 

16 

10 6 

2EC7EDB82B2DF7 



10 7 

33374B47882B32 



10 8 

3F55150F1AB3DC 



10 9 

8BA7C885CEFCE8 

log(9/10) 

10 

10 6 

80174212190900 



10 7 

21093001236414 



10 8 

01309302330968 



10 9 

44066397959215 



10 10 

82528693381274 


These computations were done at NASA Ames Research Center, using workstation 
cycles that otherwise would have been idle. 


5. Logs in base 2. 

It is easy to compute, in base 2, the d - th binary digit of 

(5.1) log(l-2-") = Ml/2"). 

So it is easy to compute log(m) for any integer m that can be written as 


(5.2) 


(2“ l - l)(2 aj - I) - •-(2 a '> - 1) 
(2*i - 1)(2 6 » - l)-..(2‘> - 1) ' 


In particular the n-th cyclotomic polynomial evaluated at 2 is so computable. A 
check shows that all primes less than 19 are of this form. The beginning of this list 


is: 

{2,3,5,7,11,13, 17,31,43,57,73,127, 151,205,257} 
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Since 


2 18 — 1 = 7 • 9 • 19 * 73, 

and since 7, y/9 and 73 are all on the above list we can compute log(19) in SC* 
from 

log(19) = log(2 18 - 1) - log(7) - log(9) - log(73). 

Note that 2 11 - 1 = 23 • 89 so either both log(23) and log(89) are in SC* or neither 
is. 

We would like to thank Carl Pomerance for showing that an identity of type (5.2) 
does not exist for 23. This is a consequence of the fact that each cyclotomic poly- 
nomial evaluated at two has a new distinct prime factor. We would also like to 
thank Robert Harley for pointing out that 29 and 37 are in SC* in base 2 via 
consideration of the Aurefeuillian factors 2 2 "" 1 + 2 n + 1 and 2 2n_1 — 2 n + C 

6. Relation Bounds. 

One of the first questions that arises in the wake of the above study is whether there 
exists a scheme of this type to compute decimal digits of 7 r. At present we know 
of no identity like (1.2) in base 10. The chances that there is such an identity are 
dimmed by some numerical results that we have obtained using the PSLQ integer 
relation algorithm [3, 12]. These computations establish (with the usual provisos 
of computer “proofs”) that there are no identities (except for the case n = 16) of 
the form 

i 00 i r 

^ _ CL\ a 2 + d 3 q m+l 

a 0 ao ^ n k mk + 1 mk -f 2 mk + m 

Jc=0 

where n ranges from 2 to 128, where m ranges from 1 to min(n,32), and where 
the Euclidean norm of the integer vector (ao,ai, • * • ,a m +i) is 10 12 or less. These 
results of course do not have any bearing on the possibility that there is a formula 
not of this form which permits computation of 7 r in some non-binary base. 

In fact, J. P. Buhler has reported a proof that any identity for 7r of the above form 
must have n = 2 A or n = \/2 . This also does not exclude more complicated 
formulae for the computation of n base 10. 

7. Questions. 

As mentioned in the previous section, we cannot at present compute decimal digits 
of 7 r by our methods because we know of no identity like (1.2) in base 10. But it 
seems unlikely that it is fundamentally impossible to do so. This raises the following 
obvious problem: 

1] Find an algorithm for the n-th decimal digit of n in SC*. It is not even clear 
that 7 T is in SC in base 10 but it ought to be possible to show this. 

2] Show that 7r is in SC in all bases. 

3] Are e and \/2 in SC (SC*) in any base? 
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Similarly the treatment of log is incomplete: 

4] Is log(2) in SC* in base 10? 

5] Is log(23) in SC* in base 2? 
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